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I 


ABSTRACT 


The  development  and  incorporation  of  the  latest  enhancements  to  the  AVAST  code  are 
described.  The  purpose  of  this  work  was  to  make  the  modelling  of  the  elasto-acoustic  interaction  more 
realistic,  while  ensuring  that  the  code  runs  as  efficiently  as  possible.  To  this  end,  several  new  features 
have  been  added.  These  include  the  option  to  use  the  wet  modes  of  the  structure  when  generating  the 
mobility  matrix,  the  option  to  generate  a  structural  impedance  matrix  based  upon  the  mass,  damping 
and  stiffness  matrices,  the  implementation  of  a  multi-frequency  boundary  integral  equation  method, 
decomposition  of  the  current  AVAST  program  modules  in  order  to  allow  for  restart  analyses' 
implementation  of  a  convergence  acceleration  technique  for  use  in  the  calculation  of  the  waveguide 
Green  s  function,  incorporation  of  the  hypersingular  Burton  and  Millar  boundary  integral  equation 
method  and  the  implementation  of  a  coupling  between  the  fluid  and  structural  models  for  the  transient 
case. 


Resume 


Cette  etude,  qui  decrit  le  developpement  et  Tincorporation  des  derniers  perfectionnements 
du  programme  AVAST,  avait  pour  but  de  rendre  encore  plus  reaiiste  la  modelisation  de  Inter¬ 
action  elasto-acoustique,  tout  en  assurant  des  passages  de  programme  aussi  efficaces  que  possi¬ 
ble.  Plusieurs  caracterisques  nouvelles  ont  ete  ajoutees  a  cette  fin,  dont  l’option  d’utilisation  des 
modes  immerges  de  la  structure  pour  la  production  d’une  matrice  d ’impedance  structural  basee 
sur  les  matrices  de  masse,  d’amortissement  et  de  rigidite  et  la  mise  en  oeuvre  d’une  methode 
d’equation  integrate  aux  limites  multi-frequences,  la  decomposition  des  modules  actuels  du  pro¬ 
gramme  AVAST  pour  permettre  les  analyses  du  redemarrage,  la  mise  en  oeuvre  d’une  technique 
de  calcul  pour  convergence  accelere  de  la  fonction  de  Green  des  guides  d’ondes,  Tincorporation 
de  la  methode  d’equation  integrate  aux  limites  hypersinguliere  de  Burton  et  Millar  et  la  mise 
en  oeuvre  d  un  couplage  entre  le  modele  fluide  et  le  modele  structural  pour  le  cas  transitoire. 
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1.1 


1.  INTRODUCTION 

Phase  I,  II  and  m  of  the  Martec/DREA  collaborative  investigation  in  underwater  acoustics  has 
resulted  in  the  development  of  a  series  of  computer  programs,  collectively  named  AVAST,  for  the 
numerical  prediction  of  the  acoustic  radiation  and  scattering  from  submerged  elastic  structures 
immersed  in  either  finite  depth,  half-space  or  infinite  fluid  domains.  AVAST  combines  both  the  finite 
element  method  for  the  structure  and  the  boundary  integral  equation  technique  for  the  surrounding 
fluid.  The  finite  element  method  is  used  to  define  a  mechanical  mobility  matrix  relating  excitation 
forces  to  structural  velocities  along  the  fluid/structure  interface.  The  boundary  integral  equation 
method  is  used  to  generate  a  system  of  equations  relating  surface  velocities  to  fluid  acoustic  pressures. 

In  an  attempt  to  improve  the  modelling  of  sound  radiated  and/or  scattered  from  submerged 
structures,  several  enhancements  have  recently  been  incorporated  into  the  previously  existing  AVAST 
suite.  These  include  the  option  to  use  the  wet  modes  of  the  structure  when  generating  the  mobility 
matrix,  the  option  to  generate  a  structural  impedance  matrix  based  upon  the  mass,  damping  and 
stiffness  matrices,  the  implementation  of  a  multi-frequency  boundary  integral  equation  method, 
decomposition  of  the  current  AVAST  program  modules  in  order  to  allow  for  restart  analyses, 
implementation  of  a  convergence  acceleration  technique  for  use  in  the  calculation  of  the  waveguide 
Green’s  function,  incorporation  of  the  hypersingular  Burton  and  Millar  boundary  integral  equation 
method  and  the  implementation  of  a  coupling  between  the  fluid  and  structural  models  for  the  transient 
case. 


In  the  discussion  which  follows,  details  concerning  the  development  and  incorporation  of  these 
latest  enhancements  to  the  AVAST  suite  will  be  presented.  In  addition,  a  state  of  the  art  literature 
review  on  infinite  wave  envelope  elements  and  a  summary  of  an  investigation  into  the  sound  radiated 
from  partially  submerged,  point  excited  bodies,  will  also  be  presented. 
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2.1 


2-  RESTRUCTURING  THE  AVAST  CODF. 

2.1  Introduction 

In  the  original  version  of  the  AVAST  suite,  all  relevant  input  data,  including  both  the  geometry 
of  the  radiating  body  and  the  position  of  the  exterior  field  points,  had  to  be  supplied  by  the  user  at  the 
start  of  an  analysis.  Once  this  information  was  collected,  the  program  would  then  proceed  with  the 
generation  of  the  associated  numerical  model  of  the  coupled  fluid/structure  system,  which  included  the 
formulation  of  frequency  dependent  fluid  matrices,  the  generation  of  a  mechanical  impedance  (or 
mobility)  matrix,  the  coupling  of  the  fluid  and  the  structural  equations,  solution  of  the  wet  surface 
velocity  and  acoustic  pressure  fields  and  finally,  the  solution  of  exterior  field  pressures.  One  major 
drawback  with  a  program  structure  of  this  type  is  the  penalty  paid  by  the  user  when  pressures  at 
additional  field  locations  are  required.  In  such  an  event  the  entire  analysis  would  have  to  be  re-run, 
representing  the  unnecessary  duplication  of  recalculating  surface  pressures  and  velocities. 

In  light  of  this  shortcoming  in  the  structure  of  the  AVAST  suite,  a  new  modular  version  of  the 
code  was  developed.  With  this  new  program  format  intermediate  results  can  now  be  saved  and  reused, 
representing  potential  savings  in  analysis  time. 

2.2  Modular  Code  Structure 


The  original  AVAST  code  has  now  been  decomposed  into  seven  distinct  program  modules: 
FLUMOD,  INWAVE,  UNCOUP,  MOBIL,  COUPLE,  PVSURF  and  EXTFLD.  The  user  may  execute 
any  or  all  modules  by  responding  positively  to  a  prompt  offered  by  the  main  driver  program.  A 
description  of  the  function  and  data  requirements  for  each  program  module  is  provided  in  separate 
sections  below. 
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2.2 

2.2.1  FLUMOD  Program  Module 

The  FLUMOD  module  generates  the  H  and  G  matrices  which  represent  the  numerical  model 
of  the  acoustic  fluid  medium  [1].  In  terms  of  input  data  requirements,  the  user  must  supply  data 
related  to  the  physical  properties  of  the  fluid,  including: 

The  number  of  structural  bodies  immersed  in  the  fluid; 

Bern  file  prefixes  (the  bem  file  contains  the  geometry  of  the  structural  wet  surface); 

The  number  of  integration  points  used  when  generating  the  H  and  G  coefficient  matrices; 

The  fluid  sound  speed  and  density; 

The  forcing  frequency;  and 

The  type  of  fluid  domain  under  consideration  (i.e.,  full  space,  half  space,  or  waveguide 

domains). 

Once  all  the  above  information  has  been  entered,  the  program  proceeds  to  generate  the  H  and 
G  coefficient  matrices,  storing  this  information  on  a  binary  data  file  having  the  extension  hgm.  In 
addition,  a  second  data  file  (having  extension  geo)  is  also  generated  by  FLUMOD  module.  This  file 
contains  data  related  to  the  geometric  properties  of  the  fluid  model  and  is  used  by  several  other 
AVAST  program  modules. 

2.2.2  INWAVE  Program  Module 

The  INWAVE  module  generates  the  incident  pressure  terms  used  in  acoustic  scattering 
analyses.  If  selected  by  the  user,  INWAVE  searches  the  working  directory  for  the  appropriate  geo  file 
generated  by  a  previous  run  of  the  FLUMOD  module.  Once  the  surface  panel  geometry  is  read  in, 
INWAVE  then  searches  the  working  directory  for  a  file  having  the  extension  src  and  sharing  the  same 
file  name  prefix  as  that  used  by  the  geo  file.  If  this  file  exists,  INWAVE  opens  it  and  extracts  the 
number  (up  to  a  maximum  of  200)  and  location  of  acoustic  sources.  The  program  then  calculates  the 
incident  pressure  values  on  the  surface  of  the  structure  and  stores  these  values  on  a  binary  data  file 
having  the  extension  inw. 


I 

I 

i 

I 

i 

I 
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2.3 


2.2.3  UNCOUP  Program  Module 


The  UNCOUP  module  allows  the  user  to  solve  for  the  wet-surface  acoustic  pressures  provided 
that  the  user  supplies  either  surface  normal  displacements,  incident  pressure  values,  or  both.  Upon 
execution,  the  program  prompts  the  user  for  the  type  of  analysis  to  be  performed,  i.e.,  radiation, 
scattering  or  combined  radiation  and  scattering.  If  an  acoustic  radiation  analysis  is  to  be  performed 
(either  with  or  without  scattering),  the  program  will  search  the  working  directory  for  a  data  file  having 
the  extension  sgd.  This  file  contains  the  wet  surface  displacements  which  are  used  to  calculate  a  set 
of  surface  acoustic  pressures.  If  an  acoustic  scattering  analysis  is  to  be  performed  (either  with  or 
without  radiation)  UNCOUP  then  searches  for  the  inw  file,  (generated  during  a  previous  run  of 
INWAVE  module),  which  contains  the  incident  source  terms.  These  incident  pressure  values  are  then 
read  in  by  the  program  and  combined,  if  necessary,  with  the  radiation  terms.  Once  the  surface  acoustic 
pressures  have  been  generated,  UNCOUP  writes  these  values  to  an  output  file  having  the  extension 
prs.  These  surface  pressure  terms  can  then  be  used,  in  combination  with  the  normal  displacements  and 
incident  pressure  terms,  to  generate  exterior  field  pressures. 

2.2.4  MOBIL  Program  Module 

The  MOBIL  program  module  generates  the  structural  mobility  matrix  required  for  coupled 
fluid/structure  interaction  problems.  If  selected  by  the  user,  MOBIL  first  searches  the  working 
directory  for  the  appropriate  VAST  [2]  compatible  t48  and  t49  files,  which  store  the  structural  stiffness 
and  mass  matrices,  respectively.  If  these  files  are  available,  MOBIL  then  prompts  the  user  for 
Rayleigh  damping  coefficients  (used  to  generate  the  structural  damping  matrix).  MOBIL  then  combines 
the  structural  mass,  stiffness  and  damping  matrices  (along  with  the  driving  frequency  found  on  the  geo 
file)  to  form  the  structural  impedance  matrix,  which  is  then  inverted  in  order  to  generate  the  structural 
mobility  matrix.  The  structural  mobility  matrix  is  then  stored  on  a  binary  data  file  having  extension 
mbl. 
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2.4 

2.2.5  COUPLE  Program  Module 

The  COUPLE  module  combines  the  structural  mobility  matrix  (stored  in  the  mbl  file)  with  the 
fluid  H  and  G  matrices  (stored  in  the  hgm  file)  to  form  a  single  system  of  equations  which  relate  fluid 
acoustic  pressures  to  structural  and  fluid  loads.  Once  generated,  this  coupled  system  of  equations  is 
then  stored  on  a  binary  data  file  having  the  extension  fsm. 

2.2.6  PVSURF  Program  Module 

The  PVSURF  module  extracts  the  coupled  system  of  equations,  stored  in  the  fsm  file,  and 
solves  for  the  wet  surface  acoustic  pressures,  surface  displacements  and  surface  normal  velocities.  The 
pressures,  displacements  and  velocities  are  then  stored  on  output  data  files  having  the  extensions  prs, 
tm  and  nrm,  respectively. 


I 

I 

I 

I 

I 

I 

I 


2.2.7  EXTFLD  Program  Module 

The  EXTFLD  module  generates  exterior  field  pressures,  provided  that  surface  acoustic 
pressures  and  normal  velocities  have  been  computed  during  a  previous  run  of  the  PVSURF  module. 
Initially,  EXTFLD  searches  the  working  directory  for  the  data  files  prs  (containing  the  surface 
pressures),  nrm  (containing  the  surface  normal  velocities)  and  efp  (containing  the  number  and  location 
of  a  series  of  exterior  field  points).  If  all  three  files  are  available  EXTFLD  will  then  proceed  with  the 
calculation  of  the  exterior  field  acoustic  pressures.  Once  this  step  has  been  completed  the  EXTFLD 
module  then  writes  these  exterior  pressures  to  an  output  data  file  having  the  extension  fpr. 
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3.1 


3.  RISC-BASED  VERSION  OF  THF,  AVAST  SUTTF. 

A  recent  review  of  DREA’s  computing  requirements  has  resulted  in  the  decision  to  replace  the 
aging  VAX  6420  VMS  mainframe  cluster  with  a  number  of  RISC-based  workstations.  As  a  result, 
much  of  the  code  now  resident  on  the  mainframe  cluster,  including  both  the  VAST  and  AVAST  suites, 
must  now  be  ported  to  these  new  platforms. 

The  installation  of  the  AVAST  code  on  the  node  RAINBOW  proceeded  without  any  difficulty. 
Previous  experience  gained  while  preparing  RISC-based  versions  of  the  code  for  the  HP-720  computer 
allowed  Martec  staff  to  generate  an  executable  version  of  the  code  on  RAINBOW  relatively  quickly. 
The  source  code  was  transferred  from  the  VAX  6420  to  the  RAINBOW  via  the  FTP  file  transfer 
facility,  where  it  was  compiled  using  the  compiler  switches  -c  and  -static. 
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4.1 

4.  STRUCTURAL  IMPEDANCE  MATRIX 

4. 1  Introduction 

In  the  original  version  of  the  AVAST  suite,  the  modal  characteristics  of  the  dry  structure  were 
used  to  establish  a  mobility  matrix,  M,  relating  surface  nodal  velocities,  v,  to  applied  structural  loads, 
f: 

(v}  =  [M]  {f}  (4.1) 

Unfortunately,  for  driving  frequencies  below  the  first  natural  frequency  of  the  dry  structure,  the  modal 
mobility  approach  failed  to  yield  acceptable  results  [3].  Experience  has  shown  that  when  the  structure 
is  excited  at  frequencies  corresponding  to  the  first  few  "wet"  natural  modes  of  the  structure  (which  are 
always  lower  than  the  first  dry  mode)  the  response  generated  by  the  mobility  approach  appears  to  be 
essentially  static. 

Since  the  structural  response  at  the  first  few  natural  frequencies  is  generally  of  significant 
interest,  a  new  structural  impedance  modelling  methodology  was  recently  proposed  [4]  and  has  now 
been  implemented  in  the  latest  version  of  the  AVAST  code.  The  impedance  modelling  approach 
provides  the  user  with  the  option  of  generating  a  structural  impedance  matrix,  Z„  based  upon  the 
structural  mass  (Ms),  damping  (Cs)  and  stiffness  (Ks)  matrices. 

PJ  -  FJ  *  io[CJ  -  <0^,]  M.2) 

where  a>  represents  the  driving  frequencies  in  radians  per  second. 

In  the  following  sections,  a  review  of  the  impedance  matrix  modelling  methodology  is 
presented.  This  is  followed  by  a  numerical  example  involving  acoustic  radiation  from  a  uniformly 
driven  pulsating  sphere  submerged  in  sea  water. 
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4.2 


4.2  Modelling  Radiated  Noise  Using  the  Impedance  Matrix  Approach 

For  the  analysis  of  time  harmonic  acoustic  radiation,  AVAST  uses  VAST  [2]  as  a  basis  to 
establish  a  finite  element  model  of  the  structure.  AVAST  combines  the  assembled  structural  stiffness, 
damping  and  mass  matrices  to  form  a  structural  impedance  matrix  in  the  form  given  above  in  Equation 
(4.1).  In  this  form,  the  impedance  matrix  provides  a  relationship  between  structural  displacements, 
x,  fluid  pressures,  p,  and  applied  structural  loads,  f, 

[Zs]{xl  =  [L](pl  +  {f}  (4-3) 

where  the  matrix  [L]  provides  a  transformation  between  surface  pressure  at  the  fluid/structure  interface 
and  structural  loads  [5]. 


The  AVAST  formulation  for  the  acoustic  fluid  is  based  upon  the  Helmholtz  integral 
relationships  [6].  The  scalar  wave  equation  for  the  acoustic  medium  is  replaced  by  an  integral  equation 
over  the  wetted  surface  (T)  of  the  structure,  which  in  turn  is  approximated  by  a  system  of  algebraic 
equations. 


[A]{pl  =  po)2[B]{xa}  +  {(]>;}  (  ) 

where  p  represents  the  fluid  density,  w  is  the  angular  frequency,  {p}  and  {xj  are  column  vectors 
whose  n-th  components  are  the  pressure  and  normal  displacement  on  the  n-th  subdivision  Tn  of  T  and 
{©.}  is  a  column  vector  whose  n-th  component  is  a  function  of  the  incident  pressure  at  a  point  on  Fn. 


In  a  pure  radiation  problem  the  incident  pressure  vanishes  and  coupling  of  the  governing 
structural  and  fluid  equations  leads  to  a  partitioned  system  for  displacement  (x)  and  pressure  (p), 


-po)2BH 


(4.5) 


where  the  matrix  [H]  provides  a  transformation  between  surface  normal  and  global  displacements. 
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4.3 


Direct  solution  of  Equation  (4.4)  is,  unfortunately,  prohibitively  expensive  for  models  of  any 
significant  size.  However,  by  eliminating  the  structural  displacements,  x.  Equation  (4.4)  can  be 
reduced  to  a  single  combined  matrix  equation  of  the  form, 

([A]  +  [R][L]){p}  =  [R]{f}  (4.6) 

where 

[R]  -  -pm2[B][H][Z]_1  (4-7) 

Once  the  fluid  acoustic  pressures  have  been  computed  they  may  then  be  combined  with  the  wet  surface 
structural  displacements  (generated  via  Equation  4.2)  to  yield  acoustic  field  pressures  via  the  exterior 
Helmholtz  equation  [1,2, 5, 6], 

4.3  Numerical  Example 

Verification  of  the  structural  impedance  modelling  feature  was  carried  out  by  comparing  results 
produced  by  AVAST  with  solutions  reported  in  the  literature  [7],  One  such  case  involved  the  acoustic 
radiation  from  a  thin-walled  spherical  shell  submerged  in  water  and  driven  internally  with  a  spherically 
symmetric  (uniform),  time  harmonic  pressure  load.  Since  the  solution  is  also  spherically  symmetric, 
the  field  solution  depends  only  on  the  radial  distance  from  the  sphere  centre.  The  physical  properties 
of  both  the  sphere  and  surrounding  fluid  are  provided  in  Table  4.1. 

By  adopting  symmetry  boundary  conditions,  the  structural  model  was  limited  to  only  a  single 
octant  of  the  sphere.  A  total  of  96  four-noded  structural  shell  elements  and  192  triangular  fluid  panel 
elements  were  used  to  model  the  structural  and  fluid  domains,  respectively.  The  numerical  model  was 
used  in  a  series  of  trial  which  involved  a  number  of  drive  frequencies.  Figure  4.1  provides  a 
comparison  of  the  AVAST  and  the  corresponding  analytical  results. 


Surface  Pressure  (Pa) 
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4.5 


Elastic  Response  of  Submerged  Sphere 

Estimates  of  Surface  Acoustic  Pressure 


Analytical  Values  — f—  AVAST  Results  Percent  Error 


FIGURE  4.1:  Elastic  Response  of  Submerged  Sphere 


Percent  Error 
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5.1 

5.  CONVERGENCE  ACCELERATION  FOR  THE  CALCULATION  OF  THE  WAVEGUIDF 
GREEN’S  FUNCTION _ _ _ _ _ 

5.1  Introduction 


The  boundary  element  method  (BEM)  has  been  an  efficient  numerical  technique  for  the 
modelling  of  acoustic  radiation  and  scattering  from  bodies  of  arbitrary  shape  immersed  in  homogeneous 
free  space  or  half  space  fluid  domains  [8].  For  a  free  space  problem,  the  Green’s  function  used  in  the 
BEM  formulation  is  simply  the  point  source  solution.  For  a  half  space  problem  with  a  perfectly 
reflecting  plane,  the  Green  s  function  can  be  obtained  by  supplementing  the  free  space  form  with  an 
image  point  source  [8,10], 


When  the  acoustic  domain  is  bounded  by  two  reflecting  planes  (i.e.,  waveguide)  the  Green’s 
function  can  no  longer  be  expressed  in  closed  form  and  must  be  represented  as  the  sum  of  an  infinite 
number  of  terms.  For  simplicity,  each  term  may  be  considered  as  a  separate  reflected  image  of  full 
space  Green’s  function. 

Although  the  concept  of  this  type  of  image  solution  is  simple  and  straightforward  [8],  its 
convergence  rate  can  be  quite  slow.  As  a  result,  Dawson  [9]  has  recently  published  a  methodology 
which  is  reported  to  significantly  accelerate  the  computation  of  the  waveguide  Green’s  function  in  cases 
where  the  series  expansion  is  ''poorly"  convergent  [9]. 

In  the  sections  which  follow,  a  description  of  Dawson’s  algorithm,  which  has  recently  been 
incorporated  into  the  AVAST  code,  is  provided.  This  is  followed  by  an  illustrative  example  which  will 
demonstrate  the  utility  of  this  new  technique. 

5.2  Green’s  Function  Acceleration 

Dawson  [9]  has  shown  that  the  Green’s  function  for  the  homogeneous  waveguide  case  (constant 
fluid  sound  speed  and  density)  may  be  expressed  in  the  following  forms. 
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<*i  = 


S(u,v  ;0)  -S(u,v+;0) 


for  Free-Free  waveguide  boundaries 


_  S(u,v_;Vfe)  -S(u,v+;Va) 
°2  4tcH 


for  Free-Rigid  waveguide  boundaries 


S(u,v  jVfe)  -  S(u,v^;V2) . 


for  Rigid-Free  waveguide  boundaries 


g4  - 


S(u,v  ;0)  -S(u,v+;0) 


for  Rigid-Rigid  waveguide  boundaries 


where  u  and  v  are  dimensionless  variables  defined  as. 


v.  =  (z  ±  %')  — 
±  v  >  H 


H  is  the  waveguide  depth,  z  and  z'  represent  the  z-coordinates  (height  above  the  waveguide  floor)  of 
the  source  and  field  points,  respectively,  r  is  the  distance  between  the  field  and  source  points  measured 
in  the  x-y  plane,  and  S(u,v,e)  can  be  expressed  as  a  sum  of  images. 
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5.3 


„  -iH«(\/u2+v2)/c7c 

S(u,v;e)  =  - ■■  -  +T(u,v;e)  +  T(u,-v;e)  (5.7) 

\/u2+v2 

where  oj  represents  the  angular  frequency,  c  represents  the  fluid  sound  speed  and  T(u,v;e)  represents 
the  half  range  sums  of  the  form, 

.  Ho)  /  7  '  2 
i - vu 

«  C7C 

T(u,v;e)  =  2  (e27tkm)— e  -  (5-8> 

m=1  \ju2  +(v  -2n:m)z 

It  has  been  reported  in  the  literature  [9]  that  non-singular  part  of  the  Green’s  function  found 
in  Equations  (5.1)  through  (5.4)  is  poorly  convergent  when  the  location  of  the  field  point  is  relatively 
close  to  the  source  point.  However,  Dawson  has  found  that  the  half  range  sums,  provided  above  in 
Equation  (5.8),  are  amenable  to  convergence  acceleration  and  can  be  rewritten  as, 

0 o 

T(u,v;e)  =  S  e2m6m+ik(27Cm-v)  Ym(u,v)  (5-9) 

m=l 


where 


k  = 


Hen 


CTC 


and 


Ym(u,v)  . 


-ik(\/u2  -*-(27tm-v)2  -  (2um  -v)) 


\jvL2  +(2ttm-v)2 

The  function  Ym  (u,v)  can  also  be  expressed  in  terms  of  the  following  expansion, 

At(u,v)  A^v)  A^ftijV) 


Ym(U>V) 


m 


m" 


0(m  -*) 


m- 


(5.10) 


(5.11) 


(5.12) 


where  the  coefficients  A,,  A2  and  A3  have  the  following  forms. 
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A,  =  — — 


(5.13) 


=  (v-iku2/2) 
^  4tt2 


(5.14) 


.  _  (v2-u2(4+8ikv-k2u2)/8) 

^  8ir3 


Dawson  then  defined  the  following  expression, 


.  AjOijV)  A^v)  A3(u,v) 

Ym(u,v)  -  -  +  - —  +  ■  ■  — 

m  m2-l  m(m2-l) 


(5.15) 


(5.16) 


which  he  then  combined  with  Equation  (5.12),  generating  the  following  formula  for  T(u,v;e) 


T(u,v;e)  =  Y,X+T*(u,v;€)  +  2  (Y„  -Y^X' 


(5.17) 


where 


X  _  e  -2iti(k-€) 


(5.18) 


and  the  function  T*(u,v;e)  is  defined  by  the  sum 


T  *(u,v;e)  =  2  Y^X' 


(5.19) 


and  can  be  evaluated  analytically.  The  coefficients  Y^  (see  Equation  5.16)  admit  a  partial  fraction 
expansion  of  the  form 


P154193.PDF  [Page:  33  of  98] 


I 
I 

I  v  v  v  A1(u,v)-A3(u,v) 

■  Ym(u,v)  =  - - - 


AjClljV)  +A3(u,v) 
2(m-l) 


A^OijV)  -  A^v) 

2(m  + 1) 


5.5 


(5.20) 


and  the  evaluation  of  the  sum  T"  then  follows  from  the  elementary  series 


m+2  m 


-X-log(l-X) 


(5.21) 


2  -Xlog(l-X)  (5.22) 

m+2  (m-1) 


S  xm  -  _2£  +  !  log(l-X) 

m+2  (m+1)  2  X 


(5.23) 


5.3  Incorporation  of  Green’s  Function  Convergence  Acceleration  Into  the  AVAST  Code 

The  convergence  acceleration  techniques  proposed  by  Dawson,  and  described  in  Section  5.2, 
have  now  been  incorporated  into  an  upgraded  version  of  the  AVAST  code.  The  significant 
improvement  in  convergence  performance  can  be  judged  qualitatively  by  considering  the  data  presented 
in  Figure  5.1,  where  Green’s  function  estimates  are  plotted  against  the  corresponding  number  of  terms 
used  in  the  series  formulation.  It  is  clear  from  the  figure  that  the  accelerated  technique  converges 
much  more  rapidly  than  does  the  original  formulation  and  hence,  representing  a  significant  potential 
saving  in  computational  time. 
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6.1 


6.  FREQUENCY  INTERPOLATION 

6.1  Introduction 

Cost-effectiveness  of  numerical  schemes  is  a  primary  consideration  in  the  development  of  any 
computational  algorithm,  since  it  determines  both  the  economics  of  computation  and  the  limit  of 
problem  solving  capabilities  [11].  For  those  researchers  involved  in  the  modelling  of  time  harmonic 
acoustic  radiation  and  scattering,  both  the  finite  element  and  boundary  element  methods  have  been 
adopted  as  the  techniques  of  choice  for  handling  large-scale  problems  [12].  The  boundary  element 
method  (BEM)  has  long  been  a  very  powerful  numerical  technique  for  acoustical  analysis,  especially 
for  exterior  acoustic  radiation  and  scattering.  The  major  advantage  of  the  BEM  over  other  numerical 
methods  is  that  only  the  surface  of  the  radiating  or  scattering  body  has  to  be  discretized  (the 
Sommerfeld  boundary  condition  is  automatically  satisfied).  However,  one  potential  shortcoming  of  the 
BEM  in  acoustics  is  that,  unlike  the  stiffness  and  mass  matrices  generated  by  the  finite  element  method, 
the  BEM  coefficient  matrices  are  frequency  dependent.  Consequently,  for  each  different  frequency  a 
new  coefficient  matrix  must  be  regenerated  [11].  This  reformulation  can  be  extremely  expensive  due 
to  the  significant  numerical  integration  involved  in  setting  up  the  matrix  equations. 

In  an  attempt  to  reduce  computation  times  for  multiple  frequency  analyses,  Benthien  and 
Schenck  [11]  have  introduced  a  so-called  "frequency  interpolation  technique."  The  strength  of  the 
method  lies  in  its  ability  to  neutralize,  or  smooth  out,  the  rapidly  varying  Green’s  function  found  in 
the  integral  equations.  Linear  interpolation  of  the  neutralized  coefficient  matrices  between  two 
reasonably  spaced  "key"  frequencies  can  then  be  performed  in  the  frequency  domain. 

In  the  sections  which  follow,  a  description  of  the  Benthien-Schenck  frequency  interpolation 
algorithm,  which  has  recently  been  incorporated  into  the  AVAST  code,  is  provided.  This  is  followed 
by  an  illustrative  example  which  will  demonstrate  the  utility  of  this  new  technique. 
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6.2  Frequency  Interpolation  Algorithm 

In  the  conventional  boundary  element  method  the  governing  differential  equation  is 
reformulated  into  a  boundary  integral  equation  defined  on  the  boundary  surface  0 


4»(p)  = 


m ,  3(j)  ,  ,  5Y 
Y(q,p)-j-(q)  -  — 
on  dn 


(q,p)  4>  (q))  dQ(q) 


where  p  and  q  represent  the  spatial  location  of  the  field  and  source  points  respectively,  n  represents 
the  surface  normal  found  at  point  q  on  the  body,  4>  represents  the  acoustic  pressure,  and  'F  represents 
the  Green’s  function,  which  for  an  infinite  fluid  domain  can  be  expressed  as: 

Y(p,q)  =  ^ -  (6-2> 


where  k  represents  the  acoustic  wavenumber  and  rpq  represents  the  distance  between  points  p  and  q. 
A  numerical  solution  of  the  integral  equation  provided  in  expression  (6.1)  can  be  achieved  by 
discretizing  the  body  surface  into  a  series  of  two-dimensional  surface  panels.  As  a  result.  Equation 
(6.1)  may  be  transformed  into  a  system  of  linear  algebraic  equations. 


[A]  V  =  [B]{“ 

On 


where, 


amj03 


=  f  A 
J  dnj 


a  -ikr(m,q) 

r(m,q) 


ds(q) 


r  -ikr(m,q) 

b-(k)  ■ 
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Examination  of  Equations  (6.4)  and  (6.5)  reveals  that  both  components  amj  and  bmj  involve  the 
quantity  e-kr,  which  varies  rapidly  with  frequency  when  the  exponent  ikr  is  relatively  large.  This  rapid 
variation  in  the  integrands  of  a^  and  bmj  can  be  eliminated  by  factoring  out  the  quantity  e  -lkRmj ,  where 
Rmj  is  ^  distance  between  centroids  of  panels  m  and  j.  Thus,  a  new  set  of  matrices  A  and  B  may 
be  defined  as 


a 


-ik[r(m,q)  -  Rmj] 

r(m,q) 


ds(q) 


(6.6) 


/ 


e  -ik[r(m,q)  -  Rmj] 

r(m,q) 


ds(q) 


(6.7) 


The  matrix  coefficients  defined  above  are  slowly  varying  with  frequency  and  are  related  to  matrices 
A  and  B  by 


%  =  e“ikR-amj(k)  (6-8) 

bmj  -e-ikR”ibmj(k)  (6-9) 

Linear  interpolation  of  the  neutralized  coefficient  matrices  can  then  be  performed  between  two 
reasonably  spaced  "key"  wavenumbers  (i.e.  kt  and  kj)  to  yield  the  matrices  A  and  B  at  each 
intermediate  wavenumber  k 


where. 


amj(k)  =  «iamj(ki)  +  MmjC^ 
hmj(k)  =  +  «2bmj(k2) 


(6.10) 

(6.11) 
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6.4 


«2  = 


(  k  -  kt  \ 

K  ~  ki, 


-ikR„ 


6.3  Incorporation  of  Frequency  Interpolation  into  the  AVAST  Code 

The  frequency  interpolation  technique  proposed  by  Benthien  and  Schenck  has  now  been 
incorporated  into  an  updated  version  of  the  AVAST  code.  In  order  to  evaluate  the  performance  of  this 
new  method  the  predicted  radiated  sound  pressure  levels  generated  by  a  uniformly  pulsating  spherical 
shell  of  radius  1.0  m,  and  computed  using  linear  interpolation,  were  compared  to  the  corresponding 
analytical  solution  found  in  the  literature.  The  acoustic  pressure  values  shown  in  Figure  6.1  were 
computed  using  base  wavenumbers  of  1.0  and  3.0  and  compare  very  favourably  with  the  exact 
solution. 


Real  Component  of  Surface  Pressure 
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FIGURE  6.1:  Acoustic  Radiation  Predicted  using  the  Frequency  Interpolation  Technique 
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7-  APPLICATION  OF  THE  BURTON  AND  MILLAR  METHOD  TN  AVAST 

7.1  Introduction 

It  has  been  well  documented  in  the  literature  that  all  classical  integral  equation  formulations 
of  the  exterior  Helmholtz  problem  fail,  either  through  non-existence  or  non-uniqueness  of  their 
solutions  at  certain  values  of  the  acoustic  wavenumber  k.  A  comprehensive  survey  of  the  most 
significant  approaches  which  have  been  proposed  to  overcome  these  mathematical  difficulties  was 
carried  out  recently  by  Burton  [13].  In  this  review  it  was  concluded  that  the  most  suitable  formulation, 
valid  for  all  wavenumbers,  is  the  one  originally  developed  by  Burton  and  Millar  [14].  The  justification 
for  this  recommendation  was  based  upon  three  principal  criteria:  reliability,  numerical  stability,  and 
efficiency.  A  second  method  proposed  by  Schenck  [15]  and  commonly  referred  to  as  the  CHIEF 
method  (Combined  Helmholtz  Integral  Equation  Formulation),  is  by  far  the  most  computationally 
efficient  scheme  available  for  exterior  acoustic  analysis,  however  it  does  not  completely  satisfy  the 
reliability  criterion  listed  above.  A  second  major  drawback  of  the  CHIEF  method  is  that  it  leads  to 
a  overdetermined  system  of  equations.  As  a  result,  for  large  scale  problems,  use  of  the  CHIEF  method 
may  be  impractical. 

Because  of  the  above  mentioned  benefits  of  the  Burton  and  Millar  method,  Martec  Limited  was 
tasked  to  incorporate  it  as  an  alternate  to  the  method  of  Schenck  in  an  upgraded  version  of  the  AVAST 
code.  In  the  sections  which  follow,  a  brief  description  of  the  method  is  provided.  This  is  followed 
by  an  illustrative  example  which  will  demonstrate  the  utility  of  this  technique. 

7.2  The  Burton  and  Millar  Integral  Formulation 

The  method  of  Burton  and  Millar  is  based  upon  a  linear  combination  of  the  surface  Helmholtz 
equation, 
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7.2 


dS(q) 


(7.1) 


and  its  normal  derivative  form 


1— <Wp)  =  f  {(Mq)— ^-(G(p,q»  -  — — -G(p,q)— ^-(|)(q)i  dS(q)  P-2> 
2  cnp  J  |  3np3aq  3np  3nq  j 


where, 

p  and  q  represent  the  acoustic  field  and  source  points  respectively; 

4>(p)  represents  the  acoustic  pressure  at  point  p; 

S  represents  the  surface  of  the  radiating  body; 

represents  the  surface  normal  at  point  q  of  the  body; 
rip  represents  the  surface  normal  at  point  p  of  the  body;  and 

G(p,q)  represents  the  Green’s  functions  appropriate  for  the  fluid  domain  under  consideration. 
This  combination  leads  to  the  following  system  of  algebraic  equations: 

([A]  *  me]  -  =  ([B]  +  pod]  -  f  (73) 


where. 
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Cij 


G(i , j)  dS.s 


dS. 


m 

13 

Sj 

ni 


represents 

represents 

represents 

represents 

represents 

represents 


the  identity  matrix 

a  coupling  parameters  (Burton  and  Millar  suggest  setting  (3  =  1/k,  where  k 
the  wavenumber) 

the  surface  of  boundary  element  panel  j 
the  surface  normal  defined  in  panel  j 
the  surface  normal  defined  in  panel  i 


While  the  uniqueness  of  the  solutions  generated  by  the  Burton  and  Millar  method  for  all 
wavenumbers  is  ensured,  a  highly  singular  kernel  (due  to  the  second  derivative  of  the  Green’s  function) 
makes  numerical  implementation  difficult.  This,  according  to  Reut  [16],  was  a  principal  reason  for 
the  method  never  gaining  wide  popularity.  Fortunately,  the  following  analytical  expressions  (exact  for 
a  flat  rectangle  or  triangle  surface  panel  Sj  containing  a  collocation  point  Xj)  were  published  in  a  recent 
article  by  Terai  [17], 


(7.4) 


and 


* 

« 


s 


j 


a2 

an..3n. 


— - —  f  — e~ikpd0, 
4 it  J  p.  J 

o  J 


(7.5) 


where  p-;  represents  the  distance  from  the  centroid  of  panel  j  to  the  panel  edge.  By  adopting  these 
definitions,  the  highly  singular  surface  integral  is  converted  into  a  non-singular  contour  integral, 
making  the  overall  numerical  implementation  straight  forward. 
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7.4 

7.3  Numerical  Example 

The  method  of  Burton  and  Millar,  incorporating  the  analytical  results  published  by  Terai,  have 
been  implemented  numerically  into  the  AVAST  code.  Numerical  tests  have  shown  that  for  a  given 
surface  mesh  the  best  accuracy  is  obtained  by  setting  the  coupling  parameter  /3  equal  to  i/k.  In  terms 
of  an  example.  Figure  7.1  provides  a  comparison  between  the  Burton  and  Millar  technique,  the 
conventional  boundary  integral  equation  formulation,  and  the  analytical  solution  for  the  acoustic 
pressure  generated  on  the  surface  of  a  uniformly  pulsating  spherical  shell.  For  driving  frequencies 
close  to  that  of  the  characteristic  frequency  (it)  the  Burton  and  Millar  formulation,  in  contrast  to  the 
conventional  method,  is  able  to  accurately  predict  surface  acoustic  pressure  levels. 
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FIGURE  7.1:  Performance  of  the  Burton  and  Millar  Method  near  a  Characteristic  Frequency 
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8.  TRANSIENT  FLUID/STRUCTURAL  COUPLING 

8.1  Introduction 

During  Phase  HI  of  the  AVAST  suite  development  work,  a  transient  response  capability  was 
outlined  and  later  implemented  into  a  new  version  of  the  code.  The  computational  method  developed 
under  this  task  combined  a  hybrid  explicit  and  implicit  time  integration  method  based  on  the  Kirchhoff 
integral  Equation  [18].  Compatibility  of  the  displacements  (defined  by  the  user)  along  the  fluid- 
structure  interface  was  guaranteed  by  the  application  of  the  Neumann  boundary  condition. 

Since  prescribing  the  displacements  is  not  practical  in  general,  the  transient  acoustic  fluid  model 
described  above  has  now  been  coupled  with  a  dynamic  model  of  an  elastic  structure.  In  the  sections 
which  follow,  a  review  of  the  fluid  hybrid  integration  method  is  provided.  This  is  followed  by  a 
description  of  the  fluid/structure  coupling  algorithm  and  the  presentation  of  an  illustrative  example. 

8.2  The  Hybrid  Time  Integration  Method 

The  acoustic  radiation  from  a  vibrating  body  is  governed  by  the  wave  equation.  The  general 
form  of  the  inhomogeneous  scalar  wave  equation  [19]  for  the  pressure  p  at  a  point  (r,t)  is  given  by, 

Ap(r,t)  -  A-^p(r,t)  =  -F(r,t)  (8.1) 

c2  at2 

with  a  Newman  type  boundary  condition  relating  normal  structural  displacements  (w)  to  fluid  pressures. 

|£(r,t)  =  ~9~  w(r,t)  (8.2) 

an  dt2 

where  F  describes  a  general  acoustic  source,  n  is  an  outward  normal  to  the  structural  boundary,  p 
represents  the  fluid  density,  c  represents  the  fluid  sound  speed  and  A  represents  the  Laplacian  operator. 
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By  using  the  time-retarded  Green’s  function  [19],  G=5[R/c the  pressure  field 
p(r,t)  may  be  described  by  the  following  integral  equation, 

c-Pfc,) .  _5_[r  1®  dS  ♦  f  -L*pt  ds 

0  K  2ir  J  s  rI  3n  I  J  s  R2  an  ^ 


r  ±2L(£)  ds  f 
s  cr  an  (  0t  Jt  J  a  R 


.If 

r2Jsl  dt  dt 


where, 

=  cosy,  —  =  _pw,  tret  =  t-R/c  (8-4) 

6n  an 

where  tQ  is  the  initial  time  and  R  is  the  distance  between  the  field  and  source  points  with  C0 — 2  for  field 
points  outside  the  body,  C0=l  for  field  points  on  the  surface  of  the  body,  and  Q=0  for  field  points 
inside  the  body.  The  symbol  Q  represents  the  exterior  fluid  domain. 

8.3  Incorporation  into  AVAST 


In  order  to  solve  Equation  (8.3),  the  surface  S  must  be  discretized  into  N  three-noded  planar 
elements  small  enough  to  ensure  that  the  variable  p  and  3p/3n  can  be  assumed  to  be  spatially  constant 
on  each  element  but  varying  with  time.  Fortunately,  models  which  have  been  developed  by  the 
AVAST  pre-processors  (stored  on  files  having  the  extension  bem)  can  be  used  by  the  new  transient 
acoustic  analysis  code. 
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Discretization  of  the  structural  surface  leads  to  the  following  form  of  Equation  (8.3), 


N 


coP(M*)  =  2  -i- 

j=l  iTZ 


r^-f-pW;)  dS.  +  f  —  cosY.-.P,  dS; 

JSjR..V  J/t"  j  J  Sj  R 2  y  tret  J 


r  1 

*1 


dS: 


*  -L  f  fM^)dq 


(8.5) 


f 

fo®-piG'| 

J  Q0 

,  *o  ftoj 

2tc  j  R 

dQ0 

!*0  =0 

with  i  describing  the  field  point,  j  the  source  point  and  R;.  =  |r.  -r.| .  If  the  source  term,  F,  in 
Equation  (8.5)  can  be  represented  by  a  Dirac  delta  function,  i.e., 

F(r,t)  =  f(t)5(r-rs)  (8.6) 

where  Rs  is  the  location  of  the  source,  then  the  integral  source  term  in  (8.5)  may  be  rewritten  as: 


L J  Ffa-RIc)  =  J_ 


R 


2  it  R_ 


(8.7) 


By  assuming  that  the  body  is  initially  at  rest,  the  last  integral  term  in  Equation  (8.5)  will 
become  zero  by  virtue  of  causality  [17],  Hence  the  pressure  fluid  at  element  i  on  the  surface  of  the 
body  can  be  expressed  at  the  k-th  time  step  in  the  form  given  below: 


N 


+  ^COSY'*~ 


+  db0SY« 


R„ 

fe] 

v  dt  It 


2tz  R_ 


(8.8) 
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8.3.1  Discretization  of  the  Time  Domain 


Assuming  that  the  time,  tret,  exists  between  tk  and  tk_l7  then  the  pressure  and  its  time  derivative 
on  element  j  may  be  approximated  in  the  form, 


“A 


k-l 


(8.9) 


I 

I 

I 

I 

I 


5 1  =  pik-pr  (8-io> 

at  ^  At 


where  a..  =  (^“^/(tk'Vi)-  Thus  whenever  At  ;>  Ry/c ,  the  pressure  at  a  field  point  on  the  j-th 
element  is  influenced  by  the  pressure  on  element  i  at  the  k-th  time  step.  This  can  be  expressed 
mathematically  as  follows: 


P(ri4)  = 


2tz 


+  ^C0SYij  r  °^Pjk 

Ky 


R, 


k-I 


^pk_pk-lN 

rj  rj 


- cosy..  , 

cRy  1J  V  At 


s.,-L5y 

J  2tt  R 


(8.11) 


8.4  Fluid/Structure  Coupling 

Compatibility  of  displacements  along  the  fluid/structure  interface  is  guaranteed  by  the 
application  of  the  Neumann  boundary  condition,  given  in  Equation  (8.2),  where  the  displacements  u(t) 
can  be  found  by  solving  the  structural  dynamic  equation, 

[M]{ii(t)}  +  [C] (u(t) }  *  [K] (u(t) }  =  {f(t)}  (8'12) 


where  [M],  [C],  and  [K]  represent  the  structural  mass,  damping,  and  stiffness  matrices  respectively, 

and  (u(t)l,  (u(t)},  and  (u(t)}  represent  generalized  structural  displacements,  velocities  and  accelerations, 

and  {f(t)l  represents  the  applied  load.  In  order  to  efficiently  treat  the  structural  dynamics  equation 
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in  the  time  domain.  Equation  (8.12)  can  be  approximated  numerically  by  adopting  the  Newmark-Beta 
method  [20],  Doing  so,  yields  an  expression  of  the  form. 


where 


[Ke]{u}t  +  At  =  {R}t+At 

[KJ  =  [K]  +  a0[M]  +  aJCJ 

{K/+At  =  [L]  {p}t+At  +  (a0{u}1  +  32  {til1  +  a3{ii}t)[M] 

+  (a^u)1  +  a4(u}t  +  a5{u}1)  [C] 

{xi}t+At  =  a0({u]t+At  -  to1)  -  a2{u)t  -  a3{ii}t 
{u}t+At  =  {u}1  +  a6{u}1  +  a7{u)t+At 

[L]  represents  a  transformation  between  fluid  pressure  and  nodal  forces. 
1 


(8.13) 


ao  = 


at  = 


®2  = 


a  At2 
5 

a  At 
1 

a  At 


a3  =  —  -  1 
2a 

6  i 

a4  =  -  -  1 
a 


as  = 


At  /  5 


-  2 


2  \  a 
a6  =  At(l  -  6) 
SLj  =  5  At 


At  represents  the  time  step  interval 

a  =  0.25  (recommended  in  [20]) 

5  =  0.5  (recommended  in  [20]) 
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The  step-by-step  solution  of  the  structural  dynamics  equation  using  the  Newmark-Beta 
integration  method  proceeds  in  three  distinct  phases.  Initially,  the  vectors  {ul ,  {u) ,  and  {til  are  set 
equal  to  the  initial  conditions  specified  by  the  user.  Following  this,  the  user  must  then  supply  the  time 
step  size  (At).  In  this  second  stage  the  program  must  form  the  effective  stiffness  matrix  ([KJ)  and 
subsequently  decompose  it  using  LU  decomposition.  In  the  final  phase  of  the  analysis  the  algorithm 
generates  the  effective  loading  on  the  structure  (which  includes  the  surface  acoustic  pressures)  at  each 
time  step  and  uses  these  to  calculate  an  updated  set  of  structural  displacements,  velocities  and 
accelerations.  The  newly  generated  structural  accelerations  are  then  substituted  into  the  Newmann 
boundary  condition  relationship  (Equation  (8.2))  and  used  to  form  an  updated  set  of  surface  acoustic 
pressures. 

8.5  Transient  Acoustic  Scattering  Example 

In  order  to  demonstrate  that  the  numerical  methods  found  in  AVAST  yield  accurate  results  for 
problems  involving  transient  acoustic  scattering,  a  model  of  a  steel  cylindrical  shell  (see  Figure  8.1) 
immersed  in  sea  water  and  subjected  to  an  exponentially  decaying  explosive  charge,  was  analyzed  by 
AVAST  and  compared  to  a  similar  model  analyzed  using  the  USA  code  (see  Table  8.1  for  the  model 
properties).  While  the  AVAST  suite  uses  the  more  rigorous  "time  retarded  potential  method",  as 
opposed  to  the  "doubly  asymptotic  approximation"  employed  in  the  USA  code,  very  close  agreement 
can  be  found  when  comparing  the  results  provided  in  Figure  8.2. 
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TABLE  8.1:  Model  Properties 

Cylinder  length 

5  m 

Cylinder  radius 

1  m 

Shell  thickness 

0.01  m 

Young’s  modulus 

2.07xl0n  Pa 

Poisson’s  ratio 

0.3 

Shell  density 

7669  kg/m3 

Fluid  sound  speed 

1500  m/s 

Fluid  density 

1000  kg/m3 

Cylinder  centroid 

x  =  0.0  m,  y  =  0.0  m,  z  =  0.0  m 

Source  location 

x  =  1000  m,  y  =  0.0  m,  z  =  0.0  m 

Source  time  profile 

e-t/0.0019 
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FIGURE  8.1:  Discretization  of  Cylindrical  Shell  Used  in  the  Transient  Acoustic  Numerical  Trials 
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9.  EXPERIMENTAL/NUMERICAL  EVALUATION  OF  PARTIALLY  SUBMERGED 
CYLINDERS _ 


9.1  Introduction 

Defence  Research  Establishment  Atlantic  (DREA)  has  recently  conducted  experiments  to 
measure  the  radiated  noise  from  a  ring-stiffened  cylinder  subjected  to  a  point-driven  harmonic  load 
[21].  These  tests,  which  were  carried  out  at  DREA  Calibration  Barge  in  Bedford  Basin  during  the 
spring  of  1993,  have  been  described  in  full  in  previous  DREA  publications  [22,23], 

In  the  report  which  follows,  selected  data  collected  during  these  experimental  trials  is  compared 
to  results  generated  by  AVAST  suite  of  programs. 

9.2  Cylinder  Data 


The  ring-stiffened  cylinder  was  manufactured  at  the  Ship  Repair  Unit  machine  shops  with 
material  provided  by  DREA.  The  cylinder  consists  of  a  longitudinally  welded  steel  tube,  stiffened  with 
ring  stiffeners  of  square  cross-section  and  capped  with  thick  steel  plating.  The  endcaps  are  of  two 
pieces  with  a  central  access  hatch,  which  is  bolted  to  the  remainder  of  the  endcap  and  sealed  with  an 
O-ring.  Ring  bolts  were  welded  to  the  endcaps  at  various  positions  to  allow  handling  of  the  cylinder 
and  the  endcaps.  The  nominal  dimensions  of  the  cylinder  are  given  in  Table  9.1. 

9.3  Numerical  Modelling 

A  finite  element  model  of  the  cylinder  was  provided  by  the  Scientific  Authority  (see  Figure 
9.1).  Table  9.2  provides  a  summary  of  the  nodes  and  elements  used  to  discretize  the  structure.  The 
fluid  model,  describing  the  wet-surface  or  fluid/structure  interface,  was  generated  by  AVAST.  A  plot 
of  the  boundary  element  fluid  mesh  is  provided  in  Figure  9.2.  A  total  of  1924  three-noded,  constant 
pressure  acoustic  panels  were  used  for  this  application.  Fluid  sound  speed  and  density  were  assumed 
constant  and  were  entered  into  the  analysis  as  1500  m/s  and  1000  kg/m3,  respectively.  It  is  also 
important  to  note  that,  due  to  proximity  of  the  structure  to  the  surface  of  the  water,  the  image  method. 
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used  to  model  half-space  fluid  domains,  was  employed  in  order  to  account  for  reflections  off  the  free 
surface. 

9.4  Comparison  of  Experimental  and  Numerical  Results 

Sound  pressure  levels  were  measured  at  a  distance  of  approximately  23  m  from  the  cylinder 
centre,  (and  at  a  depth  of  12.34  m)  for  the  cylinder  floating  horizontally  with  half  of  its  volume 
submerged.  For  this  floating  trial,  the  boundary  element  mesh  used  by  AVAST  was  simply  the  lower 
half  of  the  mesh  shown  in  Figure  9.2. 

Both  AVAST  and  VAST/COUPLE/BEMAP  codes  [24,25]  were  used  to  predict  the  radiated 
noise  from  the  floating  cylinder  at  a  frequency  corresponding  to  the  measured  2-1  mode  [2].  the  sound 
pressure  level  estimates  calculated  using  these  codes  are  compared  to  experimentally  collected  data  as 
shown  in  Figure  9.3,  providing  excellent  agreement. 


TABLE  9.1:  Nominal  Dimensions  of  the  Cylinder 

Dimension 

Value 

Length 

3000  mm 

Diameter 

762  mm 

Wall  thickness 

9.5  mm 

Endcap  thickness 

76.2  mm 

Stiffener  thickness 

38.1  mm 

Stiffener  spacing 

500  mm 

TABLE  9.2:  Summary  of  Cylinder  Structural  Model 

Geometric  nodes 

962 

Displacement  nodes 

962 

Total  degrees-of-freedom 

5772 

Shell  elements 

960 

Beam  elements 

160 

FIGURE  9.1:  Exterior  Surface  of  the  Cylinder  Structural  Finite  Element  Model 
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FIGURE  9.2:  Boundary  Element  Mesh  of  Fluid  Domain 
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FIGURE  9.3:  Directivity  Pattern  for  Floating  Cylinder 
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10.  MODELLING  RADIATED  NOISE  USING  WET  MODES 

10.1  Introduction 

When  studying  the  dynamic  response  of  structures  submerged  in  a  relatively  dense  fluid,  the 
in-vacuo  structural  modes  are  often  used  as  expansion  functions  for  obtaining  the  steady  state  far-field 
pressures  in  the  frequency  domain  [26].  Unfortunately,  for  driving  frequencies  below  the  first  natural 
frequency  of  the  dry  structure,  the  in-vacuo  modal  approach  may  fail  to  yield  acceptable  results  [3], 
Experience  has  shown  that  when  the  structure  is  excited  at  frequencies  corresponding  to  the  first  few 
"wet"  natural  modes  of  the  coupled  fluid/structure  system  (which  may  be  significantly  lower  than  even 
the  first  dry  mode)  the  response  generated  by  the  in-vacuo  modal  method  tends  to  be  essentially  static. 

As  a  result,  two  separate  algorithms  appropriate  for  determining  the  acoustic  radiation  response 
using  the  wet  modes  of  the  coupled  fluid/structure  system  have  now  been  incorporated  into  the  latest 
version  of  the  AVAST  suite.  The  first  approach  is  based  on  a  report  published  by  Ettouney  [26]  and 
uses  expansion  vectors  which  are  the  vector  components  of  the  spatial  part  of  free  harmonic  vibration 
at  the  forcing  frequency  of  the  fluid-structure  system.  In  the  second  approach  the  fluid  surrounding 
the  structure  is  assumed  to  behave  as  a  hydrodynamic  fluid  added  mass. 

In  the  sections  which  follow,  a  description  of  both  wet  mode  approaches  is  provided.  This  is 
followed  by  an  illustrative  example  which  will  demonstrate  the  relative  effectiveness  of  the  wet  mode 
techniques  at  frequencies  close  to  submerged  natural  frequencies. 

10.2  Ettounev’s  Wet  Mode  Approach 

The  governing  equation  of  the  coupled  fluid/structure  system  may  be  written  as, 

FJ{»>  ‘  F4  +  F4  <10-d 

where, 

ps]  represents  the  complex  (frequency  dependent)  dynamic  impedance  matrix  of  the  structure  in- 
vacuo 
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{8}  represents  the  structural  displacements  normal  to  the  wet  surface 
{Fej  denotes  the  structural  force  vector 

|Ff|  represents  the  force  vector  due  to  fluid/structure  interaction 

By  setting  |Fe|  equal  to  zero  and  noting  that  the  fluid/structure  interaction  force  can  be 
obtained  from  the  wet-surface  fluid  pressure  {P}  by  means  of 

{Ff}  =  [A]  {P}  (10-2) 


in  which  [A]  represents  the  diagonal  wet  surface  area  matrix,  the  following  equation  of  free  vibration 
may  be  defined 


[Ks]{5}  =  [A]{P} 


(10.3) 


In  addition,  the  fluid  pressure  {P}  may  be  expressed  as  a  series  of  m  pressure  wave  mode  shapes 


{P}  =  2  {d>k} 


(10.4) 


which  in  turn  are  related  to  a  corresponding  set  of  displacement  modes  through  the  components  aK  of 
a  diagonal  modal  mass  matrix 

{^k}  =  (10.5) 

where  a>  represents  the  system  driving  frequency  and  j£Kj-  represents  the  k-th  displacement  mode 


shape. 


Substitution  of  Equation  (10.5)  into  Equation  (10.4)  yields  the  following  expression  for  surface 


pressure 


{P}  =  -q2  2  aK  {SK} 

k=0 


(10.6) 


which,  when  inserted  into  Equation  (10.3),  yields  the  symmetric  eigenvalue  problem 
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pq  -  co2cck[A]){5k}  =  0  (10.7) 

from  which  m  eigenvectors  |5K|  (the  wet  modes)  and  eigenvalues  <xK  may  be  obtained.  It  is 

important  to  note,  however,  that  these  eigenvectors  and  eigenvalues  are  unique  to  a  specific  driving 
frequency  and  as  a  result,  a  eigenvalue  problem  must  be  solved  for  each  frequency  under  consideration. 
In  addition,  numerical  trials  have  shown  that  the  Ettouney  method  (as  implemented  in  AVAST)  is  prone 
to  convergence  problems  and  as  a  result,  has  not  been  found  to  be  very  reliable. 

10.3  Fluid  Added  Mass  Wet  Mode  Approach 

In  order  to  provide  an  alternative  to  Ettouney ’s  method,  wet  modes  computed  using  the  VAST 
code  [2]  may  also  be  used  to  model  the  acoustical  response  of  submerged  structures.  The 
computational  techniques  employed  by  VAST  to  generate  these  modes  are  based  upon  potential  theory, 
i.e.  an  inviscid,  incompressible  fluid.  In  addition,  it  is  assumed  that  the  dominant  frequency 
components  characterizing  the  motion  of  the  structure  are  low-frequency  in  nature  [27]  allowing  the 
fluid  effect  to  be  modelled  as  a  hydrodynamic  added  mass.  This  condition  refers  to  cases  in  which 

^st  >  ^ac »  Where  \st  is  a  characteristic  structural  wavelength  for  the  motion  of  the  structural  surface, 
^■ac  =  c/f  is  a  characteristic  acoustic  wavelength  for  that  motion,  f  represents  the  frequency  of  motion 

and  c  represents  the  speed  of  sound  in  the  fluid. 

Once  these  wet  modes  have  been  computed  by  VAST  they  can  be  used  to  generate  the 
structural  displacements  on  the  wet  surface.  These,  in  turn,  may  then  be  used  to  calculate  the  sound 
pressure  levels  in  the  fluid. 

10.4  Numerical  Example 

Verification  of  the  wet-mode  modelling  facilities  was  carried  out  by  comparing  results  produced 
by  the  VAST  added  mass  approach  with  output  generated  by  the  conventional  AVAST  algorithm  for 
a  case  involving  the  acoustic  radiation  produced  by  a  point  driven  cylindrical  arch  (see  Figure  10.1). 
Unfortunately,  no  meaningful  results  could  be  produced  using  the  Ettouney  algorithm.  Boundary 
conditions  were  applied  to  the  model  in  order  to  restrict  the  motion  of  the  structure  to  the  y-z  plane 
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(see  Figure  10. 1).  The  physical  properties  of  both  the  structure  and  the  surrounding  fluid  are  provided 
in  Table  10.1. 

The  First  12  wet  modes  calculated  by  VAST  were  assembled  for  this  study,  the  first  three  of 
these  are  provided  in  Figures  10.2  through  10.4.  Figure  10.5  provides  an  illustration  of  the  loads 
applied  to  the  structure,  which  consisted  of  two  point  loads  having  a  magnitude  of  1  Newton  and 
located  at  the  peak  of  the  arch.  Field  pressures  computed  for  a  location  1  meter  directly  above  the  arch 
peak,  using  both  the  wet-mode  and  conventional  methods,  is  provided  in  Figure  10.6.  It  is  important 
to  note  that  for  frequencies  very  close  to  the  predicted  wet  natural  frequency  of  approximately  72  Hz, 
the  wet  mode  approach  appears  to  indicate  a  surge  in  sound  pressure  level,  while  the  conventional 
AVAST  approach  (described  in  Section  4  of  this  report)  predicted  a  virtually  flat  response  over  this 
frequency  band. 


TABLE  10.1:  Model  Properties 

Arch  Width 

10  m 

Arch  Height 

10  m 

Shell  Thickness 

0.01  m 

Young’s  Modulus 

2.07  x  10"  Pa 

Shell  Density 

7669  Kg/m3 

Poisson’s  Ratio 

0.3 

Fluid  Sound  Speed 

1500  m/s 

Fluid  Density 

1000  kg/m3 
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Sound  Pressure  Level  (dB  re  1  .OE-6) 
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Noise  Generated  by  Vibrating  Arch 
Wet  Mode  Method  vs  Impedance  Method 


-±-  Wet  Mode  — Impedance 


FIGURE  10:6:  Noise  Generated  by  Vibrating  Arch 
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11.  INFINITE  WAVE  ENVELOPE  ELEMENTS 

11.1  Introduction 

Over  the  past  several  years,  a  number  of  numerical  methods  have  been  proposed  for  the 
modelling  of  acoustical  wave  propagation  in  infinite  or  semi-infinite  fluid  domains  [28],  Originally 
domain  techniques,  such  as  the  finite  element  and  finite  difference  methods  [29],  gained  considerable 
favour  with  researchers  active  in  the  field.  Although  obvious  difficulties  arise  in  the  application  of 
these  methods  to  infinite  domains,  proponents  of  these  techniques  simply  truncated  the  model  at  an 
arbitrary  distance  and  applied  a  suitable  boundary  condition,  such  as  the  Sommerfeld  radiation 
condition  [30],  at  the  distant  boundary.  In  general,  however,  this  lead  to  models  of  sizes  which 
exceeded  the  limits  of  conventional  computing  resources,  and  as  a  result,  these  conventional  domain 
techniques  have  been  largely  abandoned  for  more  sophisticated  methods. 

At  the  present  time  most  researchers  have  chosen  to  formulate  the  acoustic  radiation  problem 
in  terms  of  an  integral  equation  over  the  surface  of  the  radiating  body  and  solve  it  using  the  boundary 
element  method  [3],  This  has  the  two  major  advantages  in  that  the  domain  of  interest  is  reduced  from 
the  three-dimensional  exterior  region  to  a  two-dimensional  one  limited  to  the  wet  surface  of  the 
radiating  body  and  that  the  Sommerfeld  radiation  condition  is  satisfied  exactly.  However,  according 
to  Fyfe  [30],  a  number  of  disadvantages  are  associated  with  the  boundary  element  method.  From  a 
practical  point  of  view,  one  drawback  is  that  the  formulation  yields  fully  populated,  complex, 
nonsymmetric  and  frequency  dependent  coefficient  matrices.  Another  well  documented  problem  is 
related  to  the  so-called  non-uniqueness  of  the  boundary  element  solutions  at  certain  characteristic 
frequencies  [1], 

In  light  of  the  deficiencies  with  the  boundary  element  method,  considerable  effort  has  been 
devoted  to  extend  the  scope  of  the  finite  element  method  so  that  it  is  better  able  to  model  the  acoustic 
far  field  [31-35].  In  one  approach,  a  conventional  finite  element  solution  in  the  near  field  is  matched 
to  an  infinite  element  representation  in  the  unbounded  outer  region  [35-37].  Because  the  radiation 
condition  on  the  outer  boundary  requires  the  existence  of  outgoing  waves  exclusively,  such  that  all 
acoustical  energy  is  radiated  outward,  an  appropriate  asymptotic  amplitude  decay  and  wavelike 
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variation  has  been  incorporated  within  the  infinite  element  formulation.  Although  this  method  yields 
valid  solutions  in  the  inner  region  [36]  it  does  not  provide  an  accurate  solution  for  the  outer  far  field 
region. 


A  second  approach  is  that  of  the  so-called  hybrid  finite  element/  boundary  element  formulation 
[38],  whereby  a  conventional  finite  element  scheme  in  the  near  field  is  matched  with  a  boundary 
element  model  of  the  outer  fluid  domain.  According  to  Astley  [36],  this  formulation  generally  requires 
an  additional  degree  of  approximation  in  that  the  transition  from  the  near  field  (finite  element)  to  far 
field  (boundary  element)  is  seldom  as  clearly  delineated  as  the  numerical  idealization  requires.  In 
addition,  this  method  necessitates  an  iterative  matching  of  impedances  at  the  transitional  interface. 

More  recently  a  new  method,  named  the  wave  envelope  technique  [39],  has  been  proposed  in 
the  literature.  The  wave  envelope  concept  is  similar  in  many  respects  to  the  infinite  element  method. 
Both  approaches  address  the  problem  of  modelling  wave  radiation  at  the  outer  boundary  of  a 
conventional  finite  element  mesh  surrounding  a  radiating  body  in  an  infinite  domain.  They  do  so  by 
incorporating  outwardly  travelling,  wave-like  behaviour  within  the  element  shape  functions.  The 
fundamental  difference  between  the  two,  however,  lies  in  the  weighting  functions  used  in  the  Galerkin 
[20]  process,  which  forms  the  theoretical  basis  for  both  methods.  The  case  of  the  infinite  element 
formulation,  the  weighting  and  shape  functions  are  identical,  yielding  symmetric  coefficient  matrices. 
In  the  wave  envelope  formulation  the  weighting  functions  are  complex  conjugates  of  the  shape 
functions.  This  no  longer  yields  symmetric  matrices  but  does  eliminate  the  complex  harmonic  terms 
from  the  element  integrands,  providing  expressions  for  the  acoustic  mass,  stiffness  and  damping  which 
are  much  simpler  to  evaluate. 

In  the  report  to  follow  a  review  of  the  latest  research  in  the  development  of  the  wave  envelope 
methods  will  be  presented.  Of  special  interest  is  a  particular  class  of  wave  envelope  element,  named 
the  infinite  wave  envelope  element,  which  shows  a  great  deal  of  promise  in  terms  of  complimenting 


the  current  AVAST  suite. 
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11.2  The  Wave  Envelope  Method 

The  wave  envelope  method  originated  in  the  study  of  in-duct  acoustical  propagation  [40]  but 
has  been  subsequently  applied  to  homogeneous  exterior  wave  problems  [41]  and  to  inhomogeneous 
exterior  problems  involving  mean  flow  gradients  [42]  and  temperature  variations  [36],  Recently,  in 
a  paper  published  by  Astley,  Macaulay  and  Coyette  [39],  the  wave  envelope  technique  was  extended 
in  two  fundamental  respects.  First,  a  mapping  was  included  in  the  element  definition.  This  permitted 
the  wave  envelope  region  to  be  extended  to  infinity  (hence  the  name  infinite  wave  envelope  element) 
and  also  provided  greater  flexibility  in  defining  the  interface  between  the  conventional  and  wave 
envelope  regions.  As  a  rule  of  thumb,  Astley  suggested  extending  the  inner  conventional  finite  element 
mesh  to  a  distance  equivalent  to  approximately  on  acoustical  wavelength  with  roughly  seven  to  10 
elements  per  wavelength.  In  addition,  the  authors  developed  a  modified  Galerkin  scheme  in  order  to 
maintain  the  boundedness  of  the  acoustic  stiffness,  mass  and  damping.  In  essence,  this  involved  the 
incorporation  of  an  additional  weighting  factor  (equivalent  to  1/r,  where  r  approaches  infinity)  into  the 
element  formulation.  The  results  produced  by  this  modified  wave  envelope  element  were  compared 
to  both  known  analytical  solutions  and  comparable  boundary  element  studies.  For  simple  two  and  three 
dimensional  cases  the  wave  envelope  technique  was  found  to  be  in  excellent  agreement  with  the 
analytical  solutions  and  appeared  to  out-perform  the  boundary  element  method  in  terms  of  computation 
effort.  The  authors  also  suggested  that  the  infinite  wave  envelope  element  method  is  particularly  well 
suited  to  coupled  acousto-structural  problems  with  compatible,  finite  element,  structural  representations . 

Cremers  and  Fyfe  [30]  have  enhanced  the  work  of  Astley  et  al.  [39]  by  incorporating  an  n-th 
order  polynomial  in  the  element  formulation.  This  created  the  idea  of  a  powerful  element  for 
modelling  acoustic  radiation.  By  the  use  of  Lagrangian  polynomials  in  the  element  shape  functions, 
an  arbitrary  number  of  terms  in  the  1/r  expansion  can  be  generated  for  modelling  the  amplitude  decay 
of  the  outgoing  travelling  wave.  According  to  the  authors,  the  implementation  of  this  element  allows 
for  a  flexible  choice  in  the  number  of  acoustical  degrees-of-freedom  in  the  radial  infinite  direction. 
The  most  significant  advantage,  however,  is  the  potential  for  eliminating  most,  if  not  all,  of  the 
conventional  finite  elements  modelling  the  acoustic  near  field.  This  may,  in  most  cases,  significantly 
reduce  both  modelling  and  computational  times. 
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12.  USER  SUPPLIED  SURFACE  IMPEDANCE 

The  numerical  evaluation  of  the  boundary  integral  representation  of  the  Helmholtz  equation 
leads  to  an  algebraic  system  of  equations  of  the  form  [1], 

[A]{p}  =  [B]  {v}  +  {P,}  (12.1) 

where  {p}  and  {v}  are  vectors  representing  the  pressure  and  normal  velocity  evaluated  at  points  located 
on  the  fluid/structure  boundary,  {Prj.  is  a  vector  containing  the  incident  pressure  field  values,  and  [A] 
and  [B]  are  complex,  frequency  dependent,  matrices. 

When  the  surface  impedance  is  known,  an  additional  relationship  between  the  pressure  and 
velocity  may  be  defined  as  [43], 

{p}  =  -[Z]{v}  (12-2) 

where,  for  local  impedances  [43],  [Z]  is  a  diagonal  matrix  whose  elements  relate  the  pressure  and 
corresponding  normal  velocity  within  a  particular  surface  element.  Substituting  equation  (11.2)  into 
(11.1)  then  allows  for  the  direct  solution  of  the  surface  normal  velocities  using  the  formula, 

([B]  +  [A]  [Z])  {v}  =  -{Pj}  (12.3) 

In  past  versions  of  the  AVAST  suite,  the  impedance  relationship  described  above  could  only 
be  produced  by  first  developing  a  finite  element  model  of  the  structure.  In  order  to  remove  this 
restriction,  and  allow  for  greater  flexibility  in  the  definition  of  surface  impedances,  the  AVAST  code 
has  been  upgraded  to  allow  for  user  prescribed  surface  impedances.  The  AVAST  user  now  has  the 
option  of  providing  local  surface  impedances  (one  for  each  fluid  panel)  in  an  ASCII  file  having  the 
extension  usi. 

This  data  may  now  be  extracted  by  the  AVAST  code  and  used  to  set  up  a  system  of  equations 
similar  to  the  form  provided  in  Equation  (12.3). 
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